home *** CD-ROM | disk | FTP | other *** search
/ C/C++ Users Group Library 1996 July / C-C++ Users Group Library July 1996.iso / listings / v_08_05 / 8n05061a < prev    next >
Text File  |  1990-04-17  |  4KB  |  167 lines
  1. listing 2
  2.  
  3. /* programs for computations of kelvin functions and their
  4.     derivatives
  5.  
  6. DEPENDENCIES:
  7.     header file complex.h required
  8.  
  9. */
  10.  
  11. #include "complex.h"
  12. #include <stdio.h>
  13. #define max(a,b) (((a)>(b))? (a): (b))
  14. #define abs(x) ((x)?  (x):-(x))
  15. /* test driver*/
  16. void main (argc,argv) int argc; char **argv;
  17. {
  18. int i, maxit=50;
  19. double x,ber,bei,beip,berp,kei,ker,keip,kerp;
  20. FILE *fileid;
  21. fileid=fopen("PLOT.KE","w");
  22. fprintf(fileid," 4 \n");
  23. for(i=1;i<=maxit;i++)
  24.     {x=i*.25;
  25.     ke(x,&ker,&kei,&kerp,&keip);
  26.     fprintf(fileid," %f %e %e %e %e\n",x,ker,kei,kerp,keip);
  27.     };
  28. fclose(fileid);
  29. fileid=fopen("PLOT.BE","w");
  30. fprintf(fileid," 4 \n");
  31. for(i=1;i<=maxit;i++)
  32.     {x=i*.25;
  33.     be(x,&ber,&bei,&berp,&beip);
  34.     fprintf(fileid," %f %e %e %e %e \n",x,ber,bei,berp,beip);    
  35.     };
  36. fclose(fileid);
  37. exit(0);
  38. }
  39.  
  40. /* print a complex number*/
  41. printc(x) struct complex *x;
  42. {printf(" %e %e ",(x->x),(x->y));return;
  43. }
  44.  
  45. /* complex exponential*/
  46. cexp( x,ans) struct complex *x,*ans;
  47. {
  48. double y,exp(),sin(),cos();
  49. struct complex c1,c2;
  50. y = exp ( x->x);
  51. c2.x= cos (x->y);
  52. c2.y= sin (x->y);
  53. CTREAL(c1,c2,y);
  54. CLET(*ans,c1);
  55. return;
  56. }
  57.  
  58. /* Kelvin functions and their derivatives
  59.  Rational approximations from Abramowitz & Stegun
  60.  section 9.9 pp.384-5*/
  61. ke(x,ker,kei,kerp,keip)double x,*ker,*kei,*kerp,*keip;
  62. {
  63. double y,z,al,sqrt(),log(),ber,bei,berp,beip;
  64. struct complex CC,cex,cdum,f,phi,theta;
  65. if(x<=8.)
  66.  {
  67.  y=x*x/64.;
  68.  z=y*y;
  69.  be(x,&ber,&bei,&berp,&beip);
  70.  al=-log(.5*x);
  71.  
  72.  *ker=al* ber+.7853981634* bei
  73.  +(((((((-.00002458*z+.00309699)*z-.19636347)*z+5.65539121)*z
  74.  -60.60977451)*z+171.36272133)*z-59.05819744)*z-.57721566);
  75.  
  76.  *kei=al* bei-.7853981634* ber
  77.  +(((((((.00029532*z-.02695875)*z+1.17509064)*z-21.30060904)*z
  78.  +124.2356965)*z-142.91827687)*z+6.76454936)*y);
  79.  
  80.  *kerp=al* berp- ber/x+beip*.7853981634
  81.  +x*((((((-.00001075*z+.00116137)*z-.06136358)*z+1.4138478)*z
  82.  -11.36433272)*z+21.42034017)*z-3.69113734)*y;
  83.  
  84.  *keip=al* beip- bei/x-.7853981634* berp
  85.  +x*((((((.00011997*z-.00926707)*z+.33049424)*z-4.65950823)*z
  86.  +19.41182758)*z-13.39858846)*z+.21139217);
  87.  return;
  88.     }
  89. /*else*/
  90.  
  91.  CMPLX(CC,-.7071067812,-.7071067812);
  92.  CTREAL(CC,CC,x);
  93.  CTHET(-x,&y,&z);
  94.  CMPLX(theta,y,z);
  95.  CADD(CC,theta,CC);
  96.  cexp(&CC,&cex);
  97.  y=1.253314137/sqrt(x);
  98.  CTREAL(f,cex,y);
  99.  *ker=f.x;
  100.  *kei=f.y;
  101.  CPHI(-x,&y,&z);
  102.  CMPLX(phi,y,z);/*cex now phi of AS*/
  103.  CMULT(cdum,f,phi);
  104.  *kerp=-cdum.x;
  105.  *keip=-cdum.y;
  106.  return;
  107. }
  108.  
  109. be(X,BER,BEI,BERP,BEIP)double X,*BER,*BEI,*BERP,*BEIP;
  110. {
  111. struct complex CC,cex,cdum,g,theta,phi;
  112. double Y,Z,sqrt(),FKER,FKEI,FKERP,FKEIP;
  113. static double PII=.3183098862;
  114. if(X<=8.)
  115.  {
  116.  Y=(X/8.);Y=Y*Y;
  117.  Z=Y*Y    ;
  118.  *BER=((((((-.00000901*Z+.00122552)*Z-.08349609)*Z
  119.  +2.64191397)*Z-32.36345652)*Z+113.77777774)*Z-64.)*Z+1.;
  120.  *BEI=Y*((((((.00011346*Z-.01103667)*Z+.52185615)*Z
  121.  -10.56765779)*Z+72.81777742)*Z-113.77777774)*Z+16.);
  122.  *BERP=X*((((((-.00000394*Z+.00045957)*Z-.02609253)*Z
  123.  +.66047849)*Z-6.06814810)*Z+14.22222222)*Z-4.)*Y;
  124.  *BEIP=X*((((((.00004609*Z-.00379386)*Z+.14677204)*Z
  125.  -2.31167514)*Z+11.37777772)*Z-10.66666666)*Z+.5);
  126.  return;
  127.  }
  128. /*else*/
  129.     CMPLX(CC,.7071067812,.7071067812);
  130.     CTREAL(CC,CC,X);
  131.     CTHET(X,&Y,&Z);
  132.     CMPLX(theta,Y,Z);
  133.     CADD(CC,theta,CC);
  134.     cexp(&CC,&cex);
  135.     Y=.3989422804/sqrt(X);
  136.     CTREAL(g,cex,Y);
  137.     ke(X,&FKER,&FKEI,&FKERP,&FKEIP);
  138.     *BER=g.x-PII*FKEI;
  139.     *BEI=g.y+PII*FKER;
  140.     CPHI(X,&Y,&Z);
  141.     CMPLX(phi ,Y,Z);
  142.     CMULT(cdum,phi,g);
  143.     *BERP= cdum.x-PII*FKEIP;
  144.     *BEIP= cdum.y+PII*FKERP;
  145.     return;
  146. }
  147.  
  148. CTHET(X,PARTR,PARTI)double X,*PARTI,*PARTR;
  149. {
  150. double Y;
  151.  Y=8./X;
  152.  *PARTI=((((.0000019*Y+.0000051)*Y*Y-.0000901)*Y-.0009765)
  153.  *Y-.0110485)*Y-.3926991;
  154.  *PARTR=((((.0000006*Y-.0000034)*Y-.0000252)*Y-.0000906)*Y*Y
  155.  +.0110486)*Y;
  156. return;
  157. }
  158. CPHI(X,PARTR,PARTI)double X,*PARTR,*PARTI;
  159. {double Y;
  160.   Y=8./X;
  161.   *PARTI=(((((-.0000032*Y-.0000024)*Y+.0000338)*Y+.0002452)*Y
  162.   +.0013811)*Y-.0000001)*Y+.7071068;
  163.   *PARTR=((((((.0000016*Y+.0000117)*Y+.0000346)*Y+.0000005)*Y
  164.   -.0013813)*Y -.0625001)*Y+.7071068);
  165.   return;
  166. }
  167.